There are only a few methods of generating multivariate non-normal distributions. Such distributions are strongly requested by the area of multivariable analysis and statistical modeling.
Fleishmann (1978), Tadikamalla (1980), Vale and Mavrelli (1983) proposed a method for generating multivariate non-normal distributions. Its disadvantage of this method is that it is hard to calculate the population fourth order moment matrix of the simulated random vector. Yuan and Bentler (1997) proposed a procedure of generating multivariate non-normal distributions with specified marginal skewness and kurtosis considering the fourth order moment matrix.
In this description, the method I that they called is the following: ξ1, . . . , ξm are independent random variables with E (ξj)=0, E (ξj2)=1, E (ξj3)=ζj, E (ξj4)=κj and ξ=(ξ1, . . . , ξm)′. ν is a random variable which is independent of ξ, E (ν)=0, E (ν2)=1, E (ν3)=γ, E (ν4)=β. Furthermore, T=(tij) are a non-random p rows m columns matrix of rank p such that TT′=Σ and m≧p. Then the random vector given by the following formula (1),X=νTξ  (1)generally following a non-elliptical distribution with Cov (X)=Σ. Cov (•) is the covariance matrix of vector X=(x1, . . . , xn)′. The marginal skewness and kurtosis of xi are given respectively by the following formula (2) and (3)
                              skew          ⁡                      (                          x              i                        )                          =                  γ          ⁢                                    ∑                              j                =                1                            m                        ⁢                                                  ⁢                                          t                ij                3                            ⁢                                                ζ                  j                                /                                  σ                  ii                                      3                    /                    2                                                                                                          (        2        )                                          kurt          ⁡                      (                          x              i                        )                          =                  β          ⁢                      {                                                            ∑                                      j                    =                    1                                    m                                ⁢                                                                  ⁢                                                                            t                      ij                      4                                        ⁡                                          (                                                                        κ                          j                                                -                        3                                            )                                                        /                                      σ                    ii                    2                                                              +              3                        }                                              (        3        )            
The X follow an elliptical distribution when ξ˜Nm(0,I) and Pr(v≧0)=1 (Nm is m-dimensional normal distributions and Pr denotes probability) (Fang, Kotz and Ng (1990)). Let ν=1, z1, . . . , zp be independent standard variables and ξ=|z1|−E(|z1|), ξj=zj, j=2, . . . , p then the X in the formula (1) follow the skew-normal distribution as defined by Azzalini and Valle (1996).
Advantage of this method above is that the population fourth order moment matrix is easily calculated. Let S be the sample covariance based on a sample of size n from X in the formula (1). Let vech(•) be an operator which transforms a symmetric matrix into a vector by picking the non-duplicated elements of the matrix, and s=vech(S). Then the asymptotic covariance matrix of s is given Γ/n, where Γ=var(vech(XX′)) (var denotes variance) and is given by the following formula (4),
                    Γ        =                              2            ⁢            β            ⁢                                                  ⁢                                          D                p                +                            ⁡                              (                                  Σ                  ⊗                  Σ                                )                                      ⁢                          D              p                              +                ′                                              +                                    (                              β                -                1                            )                        ⁢                          vech              ⁡                              (                Σ                )                                      ⁢                                          vech                ′                            ⁡                              (                Σ                )                                              +                      β            ⁢                                          ∑                                  j                  =                  1                                m                            ⁢                                                          ⁢                                                (                                                            κ                      j                                        -                    3                                    )                                ⁢                                  vech                  ⁡                                      (                                                                  t                        j                                            ⁢                                              t                        j                        ′                                                              )                                                  ⁢                                                      vech                    ′                                    ⁡                                      (                                                                  t                        j                                            ⁢                                              t                        j                        ′                                                              )                                                                                                          (        4        )            tj being the jth column vector of T=(tij) and  is the Kronecker product. Dp is the p-order Duplication matrix. Dp+ being the Moore-Penrose generalized inverse of Dp.
After that, they studied the theory of test statics (Yuan and Bentler (1999a, 1999b, 2000)). Yuan and Bentler recommended Fleishmann (1978) method for generating random vector. However, its distribution has not so wide skewness and kurtosis. In constrast to this, the Pearson distribution system has the following characteristics.                (1) The Pearson distribution system represent wide class of distributions with various skewness and kurtosis. The Pearson system includes some well-known distributions, Gamma, Beta, t-distribution, etc.        (2) The Pearson distribution system is characterized by its four moments. This characteristics is sufficient for applying Yuan and Bentler's method.        (3) Random numbers from the Pearson distribution system except type IV is generated by using random numbers from the Normal and Gamma distributions.        (4) Recently, we developed a practical approach of using almost all types of its distribution system including the type IV distribution which was difficult to implement (Nagahara (2002)).        
We introduce the Pearson distribution system. K. Pearson (1895) defined the Pearson distribution system by the following differential equation (5) with respect to the probability density function p,
                              -                                    p              ′                        p                          =                                            b              0                        +                                          b                1                            ⁢              x                                                          c              0                        +                                          c                1                            ⁢              x                        +                                          c                2                            ⁢                              x                2                                                                        (        5        )            
Nevertheless, the Pearson type IV distribution was difficult to implement due to the difficulty in computing the normalizing constant. In the previous paper (Nagahara (1990)), by deriving recursive formulas for computing it, the author overcome various difficulties. After that, we developed a practical approach of using almost all types of its distribution system including the type IV distribution for non-Gaussian filter (Nagahara (2002)). These studies lead us to use the Pearson distribution system practically.
A glossary of researchers and texts referred to this distribution system (Pearson (1914), Pearson and Hartrey (1954), Pearson (1963), Hoadley (1968), Elderton and Johnson (1969), Ord (1972), Parrish (1983), Johnson, Kotz and Balakrishnan (1994), Stuart and Ord (1994)).
The references are the followings:                [1] Azzalini, A. and Valle, A. D. (1996) “The multivariate skew-normal distribution”, Biometrika, 83, 715-726.        [2] Devroye, L. (1984) “A simple algorithm for generating random variates with a log-concave density”, Computing, 33, 247-257.        [3] Devroye, L. (1986). Non-Uniform Random Variate Generation, Springer, New York.        [4] Elderton, W. P. and Johnson, N. L. (1969). Systems of Frequency Curves, Cambridge University Press.        [5] Fang K. - T, Kotz, S. and Ng, K. W. (1990). Symmetric Multivariate and Related Distributions, Chapman & Hall, London.        [6] Fleishmann, A. (1978). “A method for simulating non-normal distributions”, Psychometrika, 43, 521-532.        [7] Hoadley, A. B. (1968). “Use of the Pearson densities for approximating a skew density whose left terminal and first three moments are known”, Biometrika, 55, 3, 559-563.        [8] Johnson, N. L., Kotz, S. and Balakrishnan, N. (1994). Continuous univariate distributions-1, 2nd ed. John Wiley.        [9] Magnus, J. R. and Neudecker, H. (1988). Matrix Differential Calculus with Applications in Statistics and Econometrics, Wiley, New York.        [10] McGrath E. J. and Irving, D. C. (1973) “Techniques for efficient Monte Carlo simulation: volume 2: random number generation for selected probability distributions”, Technical Report SAI-72-590-LJ, Science Applications, Inc., La Jolla, Calif.        [11] Muirhead, R. J. (1982). Aspects of Multivariate Statistical Theory, Wiley, N.Y.        [12] Nagahara, Y. (1999). “The PDF and CF of Pearson type IV distributions and the ML estimation of the parameters”, Statistics & Probability Letters, Vol. 43, 251-264.        [13] Nagahara, Y. (2002). “Non-Gaussian filter and smoother based on the Pearson distribution system”, Journal of Time Series Analysis, forthcoming.        [14] Ord, J. K. (1972). Families of Frequency Distributions, Griffin.        [15] Parrish, R. (1983). “On an integrated approach to member selection and parameter estimation for Pearson distributions”, Computational Statistics & Data Analysis, 1, 239-255.        [16] Pearson, E. S. (1963). “Some problems arising in approximating to probability distributions, using moments”, Biometrika, 50, 95-111.        [17] Pearson, K. (1895). “Memoir on skew variation in homogeneous material”, Phil. Trans. Roy. Soc., A, 186, 343-414.        [18] Pearson, K. (1914). Tables for Statisticians and Biometricians, 1., Cambridge University Press.        [19] Pearson, E. S. and Hartrey, H. O. (1954). Biometrika Tables for statisticians, 1., Cambridge University Press.        [20] Shorack, G. R. (2000). Probability for Statisticians, Springer, N.Y.        [21] Stuart, A. and Ord, J. K. (1994). Kendall's Advanced Theory of Statistics, Vol. 1, Distribution Theory, Sixth edit.. Edward Arnold.        [22] Tadikamalla, P. R. (1980) “On simulating non-normal distributions”, Psychometrika, 45, 273-279.        [23] Yuan K. -H. and Bentler P. M. (1997) “Generating multivariate distributions with specified marginal skewness and kurtosis”, in SoftStat'97-Advances in Statistical Software 6, (W. Bandilla and F. Faulbaum, Eds.), 385-391, Lucius & Lucius, Stuttgart, Germany.        [24] Yuan K. -H. and Bentler P. M. (1999a). “On normal theory and associated test statistics in covariance structure analysis under two classes of nonnormal distributions”, Statist. Sinica, 9, 831-853.        [25] Yuan K. -H. and Bentler P. M. (1999b). “On asymptotic distributions of normal theory MLE in covariance structure analysis under some nonnormal distributions”, Statistics & Probability Letters, 42, 107-113.        [26] Yuan K. -H. and Bentler P. M. (2000), “Inferences on correlation coefficients in some classes of nonnormal distributions”, Journal of Multivariate Analysis, 72, 230-248.        [27] Vale, D. and Maurelli, V. A. (1983). “Simulating Multivariate Nonnormal Distribuions”, Psychometrika, 48, 465-471.        
By the way, in aspect of not only test statistics but also statistical modeling, the methods to construct the approximate distributions to the empirical distributions by introducing the non-normal distributions are more requested. And, the stable estimation of multivariate non-normal distributions from small samples are also requested.